cat(paste(Sys.Date()))
## 2018-05-31
Load libraries
library(knitr)
library(data.table)
library(tidyverse)
library(minfi)
library(minfiData)
library(sva)
library(doParallel)
library(limma)
library(DMRcate)
library(IlluminaHumanMethylation450kanno.ilmn12.hg19)
library(IlluminaHumanMethylation450kmanifest)
library(RColorBrewer)
library(ChIPpeakAnno)
library(EnsDb.Hsapiens.v75)
library(biomaRt)
library(pheatmap)
library(ggplot2)
library(gridExtra)
## Run ComBat on a GenomicRatioSet to remove batch effects.
## In:
## - A GenomicRatioSet
## - An array with the same length as the number of experiments
## in the GenomicRatioSet, with the batch for each experiment.
## Out:
## - A new GenomicRatioSet, after correction of batch effects.
combat.on.grset <- function(GRset, batch){
pheno <- pData(GRset)
combat.model = model.matrix(~1, data=pheno)
Ms_combat=ComBat(dat=getM(GRset), batch=batch, mod=combat.model)
tmp = RatioSet(M=Ms_combat, metadata=pheno)
annotation(tmp) = annotation(MsetEx)
GRset = mapToGenome(tmp)
return(GRset)
}
# Use minfi quantile normalization, then filter probes and correct batch effects.
## Input:
## - A sample table
## - Remove probes with detection p-values above this cutoff (defualt is 0.01)
## - A table with non-specific probes
## - If probes with SNPS should be dropped (defualt is TRUE)
## - If we should run ComBat to correct for batch effects (defualt is TRUE)
## Output:
## - A GenomicRatioSet object, after quantile normalization and filtering of probes.
minfi.filter.and.normalize.by.tissue <- function(sample.table, detection.pval=0.01, drop.snps=TRUE, non.specific.probe.tab=NULL, run.combat=TRUE){
RGSet <- read.metharray.exp(targets = sample.table)
GRset.quantile <- preprocessQuantile(RGSet, fixOutliers = TRUE, removeBadSamples = TRUE, badSampleCutoff = 10.5, quantileNormalize = TRUE, stratified = TRUE, mergeManifest = FALSE, sex = NULL)
# Filter on detection p-values
if(detection.pval<1){
detP <- detectionP(RGSet)
# ensure probes are in the same order in the mSetSq and detP objects
detP <- detP[match(featureNames(GRset.quantile),rownames(detP)),]
# remove any probes that have failed in one or more samples
keep <- rowSums(detP < detection.pval) == ncol(GRset.quantile)
GRset.quantile <- GRset.quantile[keep,]
}
# Remove probes overlapping SNPs
if(drop.snps){
GRset.quantile <- dropLociWithSnps(GRset.quantile , snps=c("SBE","CpG"), maf=0)
}
# Remove non-specific probes
if(!is.null(non.specific.probe.tab)){
keep <- !(featureNames(GRset.quantile) %in% non.specific.probe.tab$TargetID)
GRset.quantile <- GRset.quantile[keep,]
}
if(run.combat){
GRset.quantile <- combat.on.grset(GRset.quantile, batch=pData(GRset.quantile)$Sentrix_ID)
}
return(GRset.quantile)
}
## Input:
## - A GenomicRatioSet object
## - Which tissue we are looking at (used as header)
## Output:
## - Density plot and density bean plot
minfi.qc.plot.by.tissue <- function(GRset, sample.sheet, tissue="", sleep.column="Sleep"){
densityPlot(getBeta(GRset), sampGroups = factor( sample.sheet[,sleep.column]), main=tissue)
op <- par(mar=c(5,7,4,2)+0.1)
densityBeanPlot(getBeta(GRset), sampGroups = factor(sample.sheet[,sleep.column]),
sampNames=sample.sheet$Sample_Name, main=tissue)
par(op)
}
## Uses the DMRcate algorithm to find differentially methyltaed regions.
## In:
## - A GenomicRatioSet object
## - A string describing which model to use
## - A string describing the contrast to consider
## - FDR cutoff on differentially methylated probes (default is 0.05)
## - Fold change cutoff on differentially methylated regions (default is 0.02)
## - P-value (Stouffer) cutoff on differentially methylated regions (default is 0.05)
## - lamda: bandwidth to gaussian kernel smoothing function (default is 1000)
## - C: Scaling factor for bandwidth
## Out:
## - A GRanges object with info on all DMRs.
dmrcate.dmrs <- function(GRset, sample.sheet, use.model, contrast, cpg.fdr=0.05, dmr.fc=0.02, dmr.fdr=0.05, lambda=1000, C=2){
designMatrix <- model.matrix(as.formula(use.model), data=sample.sheet)
contMatrix <- makeContrasts(contrasts=contrast, levels=designMatrix)
myAnnotation <- cpg.annotate(GRset, datatype = "array", analysis.type="differential", design=designMatrix, arraytype="450K", contrasts = TRUE, cont.matrix = contMatrix, coef=contrast, fdr = cpg.fdr)
## If no DMRs are found, dmrcate throws an error. We need to catch this so the programs doesn't crash.
## Instead we now return an empty GRanges object.
DMRs <- tryCatch(
dmrcate(myAnnotation, lambda=lambda, C=C),
error = function(e)
{ ## Of no DMRs found, return NULL
print(e$message)
return(GRanges())
}
)
if(!(class(DMRs) == "GRanges")){
results.ranges <- extractRanges(DMRs, genome = "hg19")
results.ranges <- results.ranges[abs(results.ranges$meanbetafc) > dmr.fc & results.ranges$Stouffer < dmr.fdr,]
return(results.ranges)
}
else{
return(DMRs)
}
}
## Uses limma algorithm to find differentially methyltaed CpGc/probes.
## In:
## - A GenomicRatioSet object
## - A string describing which model to use
## - A string describing the contrast to consider
## - The number of most significant CpGs to return
## Out:
## - A a table with the most significant CpGs, and the following columns:
## ID, CHR, pos, logFC, AveExpr, t, P.Value, adj.P.Val, B
diff.cpgs <- function(GRset, sample.sheet, use.model, contrast, nr.cpgs=20){
designMatrix <- model.matrix(as.formula(use.model), data=sample.sheet)
contMatrix <- makeContrasts(contrasts=contrast, levels=designMatrix)
## Annotate CpGs.
cpg.fdr=1
myAnnotation <- cpg.annotate(GRset, datatype = "array", analysis.type="differential", design=designMatrix, arraytype="450K", contrasts = TRUE, cont.matrix = contMatrix, coef=contrast, fdr = cpg.fdr)
myAnnotation.tab <- as.data.frame(lapply(myAnnotation,head, n=1e6))
rownames(myAnnotation.tab) <- myAnnotation.tab$ID
## Run limma on each CpG.
fit <- lmFit(getBeta(GRset), designMatrix)
fit2 <- contrasts.fit(fit, contMatrix)
fit2 <- eBayes(fit2)
## Format output
tt <- topTable(fit2,coef=1, number=nr.cpgs)
tt <- cbind(myAnnotation.tab[rownames(tt),c("CHR", "pos")], tt)
tt$logFC <- signif(tt$logFC, 2)
tt$AveExpr <- signif(tt$AveExpr, 2)
tt$P.Value <- signif(tt$P.Value, 2)
tt$adj.P.Val <- signif(tt$adj.P.Val, 2)
tt$B <- signif(tt$B, 2)
return(tt)
}
## Get probes included in a DMRs found by bumphunter
## In:
## - Which program was used to get DMRs ("bumphunter" or "dmrcate")
## - DMRs (bumphunter output, bump object or DMRcate output, GRranges object)
## - The GenomicRatioSet on whic bumphunter was run
## - An index decribing wgich DMR we are interested in
## - Direction of change (1 is more methylation, -1 is less methyltaion)
## Out:
## - An array with probe ids
get.dmr.probes <- function(method, dmrs, GRset, which.dmr, change.dir){
dmr <- NULL
if(method == "bumphunter"){
dmr.tab <- dmrs$table[dmrs$table$value * change.dir > 0,]
chr <- dmr.tab[which.dmr,"chr"]
start <- dmr.tab[which.dmr,"start"]
end <- dmr.tab[which.dmr,"end"]
dmr <- GRanges(seqnames = chr, ranges = IRanges(start = start, end = end))
}
if(method == "dmrcate"){
dmr <- dmrs[which.dmr,]
}
probeinfo <- makeGRangesFromDataFrame(minfi::getAnnotation(GRset), start.field = "pos", end.field = "pos", keep.extra.columns=TRUE)
GRanges.dmr <- probeinfo[from(GenomicRanges::findOverlaps(probeinfo, dmr, ignore.strand=TRUE)),]
return( names(GRanges.dmr) )
}
## Get mean beta values for all probes in a differentially methylated region (DMR).
## Input:
## - The differentially methyltaed regions, as a GRanges object
## - The normalized methylation data, as a GenomicRatioSet object
## - Information on all probes, as a GRanges object created by makeGRangesFromDataFrame
## - Which DMR to get probes for, an intege index.
## Out:
## - An array mean beta values for a DMR, over all samples.
get.dmr.betas <- function(dmrs, GRset, probeinfo, which.dmr ){
get.dmr.probes <- function(dmrs, GRset, which.dmr, probeinfo){
dmr <- dmrs[which.dmr,]
GRanges.dmr <- probeinfo[from(GenomicRanges::findOverlaps(probeinfo, dmr, ignore.strand=TRUE)),]
return( names(GRanges.dmr) )
}
betas <- getBeta(GRset)
dmr.probes <- get.dmr.probes(dmrs, GRset, which.dmr, probeinfo)
probe.betas <- betas[dmr.probes,]
dmr.betas <- apply(probe.betas, 2, mean)
return(dmr.betas)
}
## Get differences in beta values for all probes in a differentially methylated region (DMR).
## Input:
## - The differentially methyltaed regions, as a GRanges object
## - The normalized methylation data, as a GenomicRatioSet object. ASSUMES THAT THE SAMPLES ARE
## ORDERED AS FOLLOWS: subject1 control, subject1 treat, subject2 control, subject2 treat, ...
## - Information on all probes, as a GRanges object created by makeGRangesFromDataFrame
## - Which DMR to get probes for, an intege index.
## Out:
## - An with mean beta values differences for a DMR, for all subjects in the study.
get.dmr.betas.change <- function(dmrs, GRset, probeinfo, which.dmr ){
get.dmr.probes <- function(dmrs, GRset, which.dmr, probeinfo){
dmr <- dmrs[which.dmr,]
GRanges.dmr <- probeinfo[from(GenomicRanges::findOverlaps(probeinfo, dmr, ignore.strand=TRUE)),]
return( names(GRanges.dmr) )
}
betas <- getBeta(GRset)
dmr.probes <- get.dmr.probes(dmrs, GRset, which.dmr, probeinfo)
probe.betas <- betas[dmr.probes,]
probe.beta.per.subject <- do.call("cbind", lapply(seq(from=1, to=30, by=2), function(i){ probe.betas[,i+1] - probe.betas[,i] }))
dmr.betas.per.subject <- apply(probe.beta.per.subject, 2, mean)
return(dmr.betas.per.subject)
}
## Plots differnces in methylation levels between conditions, for a singe probe of a DMR.
## Input:
## - A GenomicRatioSet
## - An id to plot, either a name of a probe in the GRset, or a name of a DMR in the GRanges object.
## - A Granges object describing DRMs. (Optional, deafualt is NULL in which case the id is expected to be a probe id.)
## - A title to put in the plot. (Default is NULL, in which case the id is used as a title).
## - Which samples correspond to the sleep condition (default 1,3,5,..,29)
## - Which samples correspond to the wake condition (default 2,4,6,..,30)
plot.methylation.difference <- function(GRset, select.id, dmrs=NULL, main=NULL, sleep=seq(from=1, to=30, by=2), wake=seq(from=2, to=30, by=2)){
betas <- getBeta(GRset)
selected.betas <- NULL
# If dmrs data is provided, assume we want to plot average over a dmr.
if(!is.null(dmrs)){
dmr.probes <- get.dmr.probes(method="dmrcate" , dmrs=dmrs, GRset=GRset, which.dmr=select.id)
probe.data <- betas[dmr.probes,]
selected.betas <- apply(probe.data,2,mean)
}
else{
selected.betas <- betas[select.id,]
}
if(is.null(main)){
main <- select.id
}
plot.data <- data.frame(sleep=pheno.data$Sleep[c(sleep,wake)], subject=pheno.data$Subject[c(sleep,wake)], beta=c(selected.betas[sleep], selected.betas[wake]))
p <- ggplot(plot.data ) +
theme_classic() +
geom_violin(aes(x = sleep, y = beta, fill=sleep), trim=FALSE) +
geom_line(aes(x = sleep, y = beta, group = subject), size=0.2) +
ggtitle(main) +
xlab("")
return(p)
}
## Given a set of genomic regions, resturns genes with TSS close to those regions.
## Input:
## - A set of genomic regions, GRanges
## - Max distance between TSS and regions (defualt=5000bp)
## - Annotation. GRanges object with
## Output:
## - An array with (Ensembl) gene ids, of genes close to the given regions.
## - An array with rownumbers of genomic regions that are close to the given Ensembl genes.
annotate.regions.to.tss <- function(regions, maxgap=5000, annoData){
overlaps.anno <- annotatePeakInBatch(regions, AnnotationData=annoData, featureType="TSS", FeatureLocForDistance="TSS", select="all", output="nearestLocation", maxgap=maxgap)
overlaps.anno$ensg <- annoData$gene_id[match(overlaps.anno$feature, annoData$tx_id)]
close.to.tss.indx <- abs(overlaps.anno$distancetoFeature)<= maxgap
close.to.regions.ensg <- unique(overlaps.anno$ensg[close.to.tss.indx])
close.to.genes.region <- unique(overlaps.anno$peak[close.to.tss.indx])
peak.gene.tab <- data.frame(region=overlaps.anno$peak, ensg=overlaps.anno$ensg)
peak.gene.tab <- peak.gene.tab[close.to.tss.indx,] ## Only keep rows where probe is close to tss
peak.gene.tab <- peak.gene.tab[order(peak.gene.tab$region,peak.gene.tab$ensg),]
## Create list mapping ensg ids <- probe ids
peak.gene.map <- split(x=as.character(peak.gene.tab$region), f=as.character(peak.gene.tab$ensg))
return(list(close.to.regions.ensg, close.to.genes.region, peak.gene.map))
}
## Convert an array of Ensembl gene ids to Entrez ids. Genes with missing ids are removed.
## Input:
## - An array with ensembl ids
## - A table with columns ensembl_gene_id and entrezgene.
## Output: An array with Entrez gene ids
ensembl.set.to.entrez <- function(ensg.ids, gene.annot.tab){
entrez.ids <- gene.annot.tab$entrezgene[match(ensg.ids, gene.annot.tab$ensembl_gene_id)]
entrez.ids <- entrez.ids[!is.na(entrez.ids)]
return(entrez.ids)
}
## Get enriched GO and KEGG terms, given a set of differentially methylated regions (DMRs).
## Input:
## - A set of DMRs (GRanges)
## - An array with (Ensembl) ids of all genes that have probes on the 450K array in the promoters.
## - A table with columns ensembl_gene_id and entrezgene.
## - Max Fisher p-value for enriched terms.
## - Minimum nr of genes for enriched terms.
## Output:
## - A list of enriched terms for: GO/BP, GO/MF, GO/CC and KEGG.
dmr.go.kegg <- function(dmrs, close.to.probe.ensembl, gene.annot.tab, annoData, fisher.cutoff=0.01, min.genes = 3){
get.top.terms <- function(enrich.obj, fisher.cutoff, min.genes){
enrich.obj[enrich.obj$P.DE< fisher.cutoff & enrich.obj$DE >= min.genes,]
}
dmr.genes <- annotate.regions.to.tss(dmrs, maxgap=5000, annoData)[[1]]
close.to.probe.entrez <- ensembl.set.to.entrez(close.to.probe.ensembl, gene.annot.tab)
dmr.genes.entrez <- ensembl.set.to.entrez(dmr.genes, gene.annot.tab)
go.fisher <- goana(dmr.genes.entrez, universe = close.to.probe.entrez, species = "Hs")
kegg.fisher <- kegga(dmr.genes.entrez, universe = close.to.probe.entrez, species = "Hs", FDR = 0.05)
out.list <- list(
BP=get.top.terms(limma::topGO(go.fisher, number=100000, truncate.term=40, ontology = "BP"), fisher.cutoff, min.genes),
MF=get.top.terms(limma::topGO(go.fisher, number=100000, truncate.term=40, ontology = "MF"), fisher.cutoff, min.genes),
CC=get.top.terms(limma::topGO(go.fisher, number=100000, truncate.term=40, ontology = "CC"), fisher.cutoff, min.genes),
KEGG=get.top.terms(limma::topKEGG(kegg.fisher, number=100000), fisher.cutoff, min.genes)
)
return(out.list)
}
## Modified version of the function ChIPpeakAnno::assignChromosomeRegion.
## Instread of returning percentage and jacard index of regions annotated to different
## genomic features, an array with the feature for each input region is returned.
## Input:
## - the genomic regions to annoate (GRanges object)
## - upstream cutoff for promoters (default 1000)
## - downstream cutoff for promoters (default 1000)
## - Precedence of genomic features, in case a region overlaps several different fetures.
## Has to be any combination of "Exon", "Intron", "fiveUTR", "threeUTR", "Promoter",
## and "immediateDownstream".
## - A txdb with the genomic features.
## Output:
## - An array with the type of genomic feature for each input region.
## See ?assignChromosomeRegion for more info.
assignChromosomeRegion2 <- function (peaks.RD, proximal.promoter.cutoff = 1000L,
immediate.downstream.cutoff = 1000L, precedence, TxDb)
{
if (!inherits(TxDb, "TxDb"))
stop("TxDb must be an object of TxDb, \n try\n?TxDb\tto see more info.")
if (!inherits(peaks.RD, c("RangedData", "GRanges")))
stop("peaks.RD must be a GRanges object.")
if (!is.null(precedence)) {
if (!all(precedence %in% c("Exon", "Intron", "fiveUTR",
"threeUTR", "Promoter", "immediateDownstream")))
stop("precedence must be a combination of \nExons, Introns, fiveUTRs, threeUTRs, \nPromoters, immediateDownstream")
}
if (inherits(peaks.RD, "RangedData")){
peaks.RD <- as(peaks.RD, "GRanges")
}
exons <- exons(TxDb, columns = NULL)
introns <- unique(unlist(intronsByTranscript(TxDb)))
fiveUTRs <- unique(unlist(fiveUTRsByTranscript(TxDb)))
threeUTRs <- unique(unlist(threeUTRsByTranscript(TxDb)))
transcripts <- unique(transcripts(TxDb, columns = NULL))
options(warn = -1)
try({
promoters <- unique(promoters(TxDb, upstream = proximal.promoter.cutoff,
downstream = 0))
immediateDownstream <- unique(flank(transcripts,
width = immediate.downstream.cutoff, start = FALSE,
use.names = FALSE))
})
microRNAs <- tryCatch(microRNAs(TxDb), error = function(e) return(NULL))
tRNAs <- tryCatch(tRNAs(TxDb), error = function(e) return(NULL))
options(warn = 0)
annotation <- list(exons, introns, fiveUTRs, threeUTRs, promoters, immediateDownstream)
if (!is.null(microRNAs))
annotation <- c(annotation, microRNAs = microRNAs)
if (!is.null(tRNAs))
annotation <- c(annotation, tRNAs = tRNAs)
annotation <- lapply(annotation, function(.anno) {
mcols(.anno) <- NULL
.anno
})
names(annotation)[1:6] <- c("Exon", "Intron", "fiveUTR",
"threeUTR", "Promoter", "immediateDownstream")
formatSeqnames <- function(gr) {
seqlevels(gr)[grepl("^(\\d+|MT|M|X|Y)$", seqlevels(gr))] <- paste("chr",
seqlevels(gr)[grepl("^(\\d+|MT|M|X|Y)$", seqlevels(gr))],
sep = "")
seqlevels(gr)[seqlevels(gr) == "chrMT"] <- "chrM"
trim(gr)
}
peaks.RD <- formatSeqnames(peaks.RD)
peaks.RD <- unique(peaks.RD)
annotation <- lapply(annotation, formatSeqnames)
annotation <- GRangesList(annotation)
newAnno <- c(unlist(annotation))
newAnno.rd <- GenomicRanges::reduce(trim(newAnno))
Intergenic.Region <- gaps(newAnno.rd, end = seqlengths(TxDb))
Intergenic.Region <- Intergenic.Region[strand(Intergenic.Region) != "*"]
if (!all(seqlevels(peaks.RD) %in% seqlevels(newAnno))) {
warning("peaks.RD has sequence levels not in TxDb.")
sharedlevels <- intersect(seqlevels(newAnno), seqlevels(peaks.RD))
peaks.RD <- keepSeqlevels(peaks.RD, sharedlevels)
}
mcols(peaks.RD) <- NULL
annotation <- annotation[unique(c(precedence, names(annotation)))]
names(Intergenic.Region) <- NULL
annotation$Intergenic.Region <- Intergenic.Region
anno.names <- names(annotation)
ol.anno <- findOverlaps(peaks.RD, annotation)
ol.anno.splited <- split(queryHits(ol.anno), anno.names[subjectHits(ol.anno)])
ol.anno <- as.data.frame(ol.anno)
ol.anno.splited <- split(ol.anno, ol.anno[, 2])
hasAnnoHits <- do.call(rbind, ol.anno.splited[names(ol.anno.splited) != as.character(length(annotation))])
hasAnnoHits <- unique(hasAnnoHits[, 1])
ol.anno <- ol.anno[!(ol.anno[, 2] == length(annotation) & (ol.anno[, 1] %in% hasAnnoHits)), ]
ol.anno <- ol.anno[!duplicated(ol.anno[, 1]), ]
out.data <- names(annotation)[ol.anno$subjectHits]
}
## Annotate DMRs, with genomic region and (possible) promoter.
## Input:
## - DMRs, as GRanges object
## - Max distance to transcription start site (default 5000 bp)
## - Annotation data as ??
## - Annotation data as txdb
## Output:
## - The same GRanges object as input, but with two additional columns:
## "gene_tss" (if the DMR is close to any TSS) and "genomic _feature"
## ("near TSS", "exon", "intron", "Intergenic.Region" etc.)
annot.dmrs <- function(dmrs, maxgap=5000, annoData, annoData.TXDB){
## Find transcription start sites nearby the DMRs
overlaps.anno <- annotatePeakInBatch(dmrs, AnnotationData=annoData, featureType="TSS", FeatureLocForDistance="TSS", select="all", output="nearestLocation", maxgap=maxgap)
overlaps.anno$gene_name <- annoData$gene_name[match(overlaps.anno$feature, names(annoData))]
overlaps.anno <- overlaps.anno[abs(overlaps.anno$distancetoFeature)<= maxgap, ] ## Only keep nearby TSS
## For each DMR, get names of genes with nearby TSS (if any)
nearest.gene <- sapply(names(dmrs),
function(x){
peak.indx <- overlaps.anno$peak == x & base::grepl("ENST",overlaps.anno$feature)
paste(unique(overlaps.anno$gene_name[peak.indx]), collapse=", ")
})
dmrs$gene_tss <- nearest.gene
assigned.chromosome.regions <- assignChromosomeRegion2(dmrs, precedence=c("Promoter", "immediateDownstream","fiveUTR", "threeUTR", "Exon", "Intron"), proximal.promoter.cutoff=maxgap, immediate.downstream.cutoff=maxgap, TxDb=annoData.TXDB)
dmrs$genomic_feature <- assigned.chromosome.regions
## Make sure the the tss annotations from annotatePeakInBatch and the assigned chromatin regions agree.
gene.tss.overlap <- dmrs$gene_tss != ""
dmrs$genomic_feature[gene.tss.overlap] <- "near TSS"
return(dmrs)
}
## Summarize probe data by gene.
## Input:
## - The GenomicRatioSet with the probe data.
## - A list mapping gene id -> probe ids
## - Which function to use for summarizing. Default is mean.
## Output:
## - A matrix with gene as rows, samples as column, with summarized M-values of all probes annotated
## to the gene in question. E.g. if the mean function is used, mean M-values over all probes for a
## gene are returned.
mvals.per.gene <- function(GRset, probe.gene.map, use.function=mean){
filtered.map <- lapply(probe.gene.map,
function(x){
intersect(x, rownames(GRset))
})
mvals <- getM(GRset)
out.list <- lapply(names(filtered.map),
function(ensg){
probe.data <- mvals[filtered.map[[ensg]], ]
if(length(nrow(probe.data)) ==0){
return(NULL)
}
return( apply(probe.data, 2, use.function, na.rm=TRUE) )
})
names(out.list) <- names(filtered.map)
out.tab <- do.call(rbind, out.list)
return(out.tab)
}
## Export probe data as gzipped csv files with M-values.
## Input:
## - The GenomicRatioSet to export.
## - Which probes to export.
## - The name of the output file.
## Output:
## - Writes the M values from the given GenomicRatioSet to the given file.
write.m.vals <- function(GRset, use.probes, out.file){
out.data <- getM(GRset)[rownames(GRset) %in% use.probes, ]
out.data <- round(out.data,4)
out.file.gz <- gzfile(out.file)
write.csv(out.data , file=out.file.gz)
}
# setwd()
set.seed(1977) ## set seed so that methods using randomization (bumphunter and SAM) always give the same results.
#########################################
## Minfi load and pre-process data
## Read info on non-specific probes
non.specific.probe.tab <- fread('https://raw.githubusercontent.com/sirselim/illumina450k_filtering/master/48639-non-specific-probes-Illumina450k.csv', header=TRUE)
## Read table with info about all arrays
targets <- read.csv(file="cedernaes2018_dnamethylation_samplesheet.csv")
rownames(targets) <- paste(targets$Sentrix_ID, targets$Sentrix_Position, sep="_")
targets.a <- targets[targets$Tissue=="A",]
targets.m <- targets[targets$Tissue=="M",]
## Normalize. Since there are large differences between adipose and muscle, use quantile normalization on each tissue separately.
minfi.norm.a <- minfi.filter.and.normalize.by.tissue (targets.a, non.specific.probe.tab=non.specific.probe.tab)
## Warning in .getSex(CN = CN, xIndex = xIndex, yIndex = yIndex, cutoff
## = cutoff): An inconsistency was encountered while determining sex. One
## possibility is that only one sex is present. We recommend further checks,
## for example with the plotSex function.
## Standardizing Data across genes
minfi.norm.m <- minfi.filter.and.normalize.by.tissue (targets.m, non.specific.probe.tab=non.specific.probe.tab)
## Warning in .getSex(CN = CN, xIndex = xIndex, yIndex = yIndex, cutoff
## = cutoff): An inconsistency was encountered while determining sex. One
## possibility is that only one sex is present. We recommend further checks,
## for example with the plotSex function.
## Standardizing Data across genes
minfi.qc.plot.by.tissue(minfi.norm.a, sample.sheet=targets.a, tissue="Adipose (minfi)")
minfi.qc.plot.by.tissue(minfi.norm.m, sample.sheet=targets.m, tissue="Muscle (minfi)")
## Get DMRs with dmrcate.
dmrcate.minfi.dmrs.a <- dmrcate.dmrs(minfi.norm.a, sample.sheet=targets.a, use.model="~ 0 + Sleep + Subject", contrast="Sleepwake-Sleepsleep", dmr.fc=0.02)
dmrcate.minfi.dmrs.m <- dmrcate.dmrs(minfi.norm.m, sample.sheet=targets.m, use.model="~ 0 + Sleep + Subject", contrast="Sleepwake-Sleepsleep", dmr.fc=0.02)
## [1] "The FDR you specified in cpg.annotate() returned no significant\n CpGs, hence there are no DMRs. Try specifying a value of\n 'pcutoff' in dmrcate() and/or increasing 'fdr' in\n cpg.annotate()."
dmrcate.minfi.dmrs.nofc.a <- dmrcate.dmrs(minfi.norm.a, sample.sheet=targets.a, use.model="~ 0 + Sleep + Subject", contrast="Sleepwake-Sleepsleep", dmr.fc=0)
dmrcate.minfi.dmrs.nofc.m <- dmrcate.dmrs(minfi.norm.m, sample.sheet=targets.m, use.model="~ 0 + Sleep + Subject", contrast="Sleepwake-Sleepsleep", dmr.fc=0)
## [1] "The FDR you specified in cpg.annotate() returned no significant\n CpGs, hence there are no DMRs. Try specifying a value of\n 'pcutoff' in dmrcate() and/or increasing 'fdr' in\n cpg.annotate()."
There were no signficant DMRs in muscle. Therefore we also look at the the most signficant individual probes/CpGs. None of these are very signficant either (from the FDR adjusted p-value P.adj):
diff.cpgs(GRset=minfi.norm.m, sample.sheet=targets.m, use.model="~ 0 + Sleep + Subject", contrast="Sleepwake-Sleepsleep", nr.cpgs=10)
## CHR pos logFC AveExpr t P.Value adj.P.Val B
## cg02934600 chr1 107600376 0.028 0.34 6.774544 3.9e-06 1 4.0
## cg06070414 chr7 95226412 0.044 0.25 6.493057 6.6e-06 1 3.4
## cg00722081 chr10 61147668 -0.028 0.63 -6.327758 9.0e-06 1 3.1
## cg26493306 chr5 134463458 0.023 0.92 6.275672 9.9e-06 1 3.0
## cg04957146 chr17 2551652 0.021 0.87 6.164014 1.2e-05 1 2.8
## cg01051642 chr1 2517410 0.026 0.12 6.000548 1.7e-05 1 2.5
## cg09250087 chr2 161127564 0.029 0.69 5.862363 2.2e-05 1 2.2
## cg20292318 chr17 77141853 -0.025 0.82 -5.494034 4.5e-05 1 1.5
## cg10747185 chr3 176862267 0.038 0.84 5.493864 4.5e-05 1 1.5
## cg03063453 chr19 52693134 0.023 0.13 5.490893 4.5e-05 1 1.5
pval.cutoff <- 0.05
length(which(dmrcate.minfi.dmrs.a$meanbetafc > 0))
## [1] 62
length(which(dmrcate.minfi.dmrs.a$meanbetafc < 0))
## [1] 31
length(which(dmrcate.minfi.dmrs.m$meanbetafc > 0))
## [1] 0
length(which(dmrcate.minfi.dmrs.m$meanbetafc < 0))
## [1] 0
length(which(dmrcate.minfi.dmrs.nofc.a$meanbetafc > 0))
## [1] 92
length(which(dmrcate.minfi.dmrs.nofc.a$meanbetafc < 0))
## [1] 56
length(which(dmrcate.minfi.dmrs.nofc.m$meanbetafc > 0))
## [1] 0
length(which(dmrcate.minfi.dmrs.nofc.m$meanbetafc < 0))
## [1] 0
maxgap <- 5000
annoData <- toGRanges(EnsDb.Hsapiens.v75, feature="transcript")
annoData.TXDB <- makeTxDbFromGFF("/Users/orzechoj/projects/benedict2016/data/annot/Homo_sapiens.GRCh37.87.gff3.gz", format="gff3")
## Import genomic features from the file as a GRanges object ... OK
## Prepare the 'metadata' data frame ... OK
## Make the TxDb object ... OK
dmrcate.minfi.dmrs.a <- annot.dmrs(dmrcate.minfi.dmrs.a, maxgap=5000, annoData, annoData.TXDB)
dmrcate.minfi.dmrs.nofc.a <- annot.dmrs(dmrcate.minfi.dmrs.nofc.a, maxgap=5000, annoData, annoData.TXDB)
pheno.data <- targets.a
rownames(pheno.data) <- paste(pheno.data$Sentrix_ID, pheno.data$Sentrix_Position, sep="_")
probeinfo <- makeGRangesFromDataFrame(minfi::getAnnotation(minfi.norm.a), start.field = "pos", end.field = "pos", keep.extra.columns=TRUE)
## Plot all samples, DMRs with any fold change
plot.data <- t(sapply(1:length(dmrcate.minfi.dmrs.nofc.a), function(i){get.dmr.betas(dmrcate.minfi.dmrs.nofc.a, minfi.norm.a, probeinfo, i)}))
plot.data <- plot.data[,c(seq(from=1, to=ncol(plot.data), by=2),seq(from=2, to=ncol(plot.data), by=2))]
pheatmap(plot.data, cluster_rows=TRUE, cluster_cols=FALSE, show_colnames=FALSE, scale="row",
annotation_col=pheno.data[,c("Sleep","Subject")], main="Adipose DMRs (average Beta values)", fontsize_row=5)
## Which DMRs to print
plot.names <- c("CD36", "TSPAN", "GNAS", "INS", "GFI1", "AKR1CL1", "TNXB", "TRIM2", "HOXA2","FOXP2")
plot.dmrs <- names(dmrcate.minfi.dmrs.nofc.a[unlist(lapply(plot.names, function(x){grep(x,dmrcate.minfi.dmrs.nofc.a$gene_tss)})),] )
names(plot.dmrs) <- plot.names
## Make violin plots
p <- lapply(1:length(plot.dmrs),
function(i){
plot.methylation.difference(GRset=minfi.norm.a, select.id=plot.dmrs[i], dmrs=dmrcate.minfi.dmrs.nofc.a, main=names(plot.dmrs)[i])
})
ml <- marrangeGrob(p, nrow=5, ncol=2)
print(p)
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
##
## [[9]]
##
## [[10]]
## Fetch gene info.
ensembl = useEnsembl(biomart="ensembl",GRCh=37, dataset="hsapiens_gene_ensembl")
gene.annot.tab <- getBM(attributes=c('ensembl_gene_id','gene_biotype','entrezgene', 'external_gene_name'), mart = ensembl)
## Clean up data: if there are several entries with same ensembl_gene_id, use the first
gene.annot.tab %>%
group_by(ensembl_gene_id) %>%
summarise_each(funs(paste(sort(.), collapse=","))) %>%
extract(external_gene_name, "external_gene_name", "([^,]+),?.*") %>%
extract(gene_biotype, "gene_biotype", "([^,]+),?.*") %>%
extract(entrezgene, "entrezgene", "([^,]+),?.*") -> gene.annot.tab
## `summarise_each()` is deprecated.
## Use `summarise_all()`, `summarise_at()` or `summarise_if()` instead.
## To map `funs` over all variables, use `summarise_all()`
## Create genomic ranges object from info about 450K probes
ann450k = minfi::getAnnotation(IlluminaHumanMethylation450kanno.ilmn12.hg19)
cpgData <- GRanges(seqnames=Rle(ann450k$chr),
ranges=IRanges(start=ann450k$pos, end=ann450k$pos),
strand=Rle(rep("*",nrow(ann450k))))
names(cpgData) <- rownames(ann450k)
## Get all genes that have 450K probes in promoters
probe.gene.overlaps <- annotate.regions.to.tss(cpgData, maxgap=5000, annoData)
close.to.probe.ensembl <- probe.gene.overlaps[[1]]
probes.near.genes <- probe.gene.overlaps[[2]]
probe.gene.map <- probe.gene.overlaps[[3]]
Using topGO and topKEGG from limma.
## Get enriched pathways etc. in genes close to hypermethylated adipose DMRs
enriched.terms.a.no.fc <- dmr.go.kegg(dmrcate.minfi.dmrs.nofc.a, close.to.probe.ensembl, gene.annot.tab, annoData, fisher.cutoff=0.01, min.genes = 3)
enriched.terms.a.down.no.fc <- dmr.go.kegg(dmrcate.minfi.dmrs.nofc.a[dmrcate.minfi.dmrs.nofc.a$meanbetafc <0,],
close.to.probe.ensembl, gene.annot.tab, annoData, fisher.cutoff=0.01, min.genes = 3)
enriched.terms.a.up.no.fc <- dmr.go.kegg(dmrcate.minfi.dmrs.nofc.a[dmrcate.minfi.dmrs.nofc.a$meanbetafc >0,],
close.to.probe.ensembl, gene.annot.tab, annoData, fisher.cutoff=0.01, min.genes = 3)
Enrichr ia a web tool where you can upload lists of genes, and look for overlaps with gene sets from many different sources (Gene Ontology, KEGG, Wikipathays, ENCODE, GTEx). The tool can be found at http://amp.pharm.mssm.edu/Enrichr/. The following code generates gene lists that can be pasted into the web form.
write.table(dmrcate.minfi.dmrs.a[dmrcate.minfi.dmrs.a$meanbetafc > 0,]$gene_tss, row.names = F, col.names = F, quote = F)
## TNXB
## ABAT
## PLEKHB1
## PHKG1
## LGI1
## PRDM1
## LGALS3BP
## NRXN1
## C8orf86
## RIN2
## PCDH15
## SOD3
## FZD4
## VMP1
## BAIAP2
## PRLR
## DCN
## RIMKLB
## IFLTD1
## CHD3
## FOXP2
##
## TRIM2
## AKR1CL1
## ADAM5
## ASB16
## RP11-820L6.1
##
## SPG20
## LINC00917
## RP11-209A2.1
## SULT1C2
## SEMA5B
## HMGB1P13
## CD36
## SLC38A4
## TEX26-AS1
## CTC-529P8.1
##
## SMG5
## BCHE
##
## SYNPO
## TMEM229B
## STARD8
## PLCD4
## CES3
## RNF183
## SAPCD1-AS1
##
## SNAI2
## TMEM244
## COL11A2
##
## CCDC80
## C7orf62
## HOXA2
## FXYD3
## PLA2G3
## CHM
## RGS12
## ST5
write.table(dmrcate.minfi.dmrs.a[dmrcate.minfi.dmrs.a$meanbetafc < 0,]$gene_tss, row.names = F, col.names = F, quote = F)
## TSPAN32
##
## AVP
## WDFY4
## SNORD115-18
## ZNF385A
## TOP1MT
## RP11-627G18.2
## COL10A1
## ADORA2A
## UNCX
## MSX1
##
## MAD1L1
##
## CPT1A
## RHOD
##
## RP11-441F2.5
## SLC25A30
## HOXD3
## RGS12
## TRIM39
## TMC6
## MPG
## RP11-109A6.3
##
## HIC1
## GFI1
## MSX1
## CTSZ
write.table(dmrcate.minfi.dmrs.a$gene_tss, row.names = F, col.names = F, quote = F)
## TNXB
## TSPAN32
##
## ABAT
## PLEKHB1
## PHKG1
## AVP
## LGI1
## WDFY4
## PRDM1
## LGALS3BP
## NRXN1
## C8orf86
## RIN2
## PCDH15
## SOD3
## FZD4
## VMP1
## SNORD115-18
## BAIAP2
## PRLR
## DCN
## ZNF385A
## RIMKLB
## IFLTD1
## CHD3
## FOXP2
##
## TRIM2
## AKR1CL1
## ADAM5
## TOP1MT
## ASB16
## RP11-820L6.1
##
## SPG20
## RP11-627G18.2
## LINC00917
## COL10A1
## RP11-209A2.1
## SULT1C2
## SEMA5B
## ADORA2A
## HMGB1P13
## CD36
## SLC38A4
## UNCX
## TEX26-AS1
## MSX1
## CTC-529P8.1
##
##
## MAD1L1
##
## SMG5
## CPT1A
## RHOD
##
## RP11-441F2.5
## SLC25A30
## BCHE
##
## SYNPO
## HOXD3
## TMEM229B
## STARD8
## PLCD4
## CES3
## RNF183
## SAPCD1-AS1
## RGS12
## TRIM39
##
## SNAI2
## TMEM244
## TMC6
## COL11A2
##
## MPG
## CCDC80
## C7orf62
## HOXA2
## RP11-109A6.3
## FXYD3
## PLA2G3
##
## HIC1
## GFI1
## MSX1
## CTSZ
## CHM
## RGS12
## ST5
Uploading the lists to Enrichr gave the following results:
Adipose down http://amp.pharm.mssm.edu/Enrichr/enrich?dataset=3hvym
Adipose up http://amp.pharm.mssm.edu/Enrichr/enrich?dataset=3hvxp
Adipose all http://amp.pharm.mssm.edu/Enrichr/enrich?dataset=3hvz1
write.table(as.data.frame(dmrcate.minfi.dmrs.a[dmrcate.minfi.dmrs.a$meanbetafc > 0,]), sep="\t", quote=F, row.names=F, file="./minfi_dmrcate_a_up.txt")
write.table(as.data.frame(dmrcate.minfi.dmrs.a[dmrcate.minfi.dmrs.a$meanbetafc < 0,]), sep="\t", quote=F, row.names=F, file="./minfi_dmrcate_a_down.txt")
write.table(as.data.frame(dmrcate.minfi.dmrs.nofc.a[dmrcate.minfi.dmrs.nofc.a$meanbetafc > 0,]), sep="\t", quote=F, row.names=F, file="./minfi_dmrcate_nofc_cutoff_a_up.txt")
write.table(as.data.frame(dmrcate.minfi.dmrs.nofc.a[dmrcate.minfi.dmrs.nofc.a$meanbetafc < 0,]), sep="\t", quote=F, row.names=F, file="./minfi_dmrcate_nofc_cutoff_a_down.txt")
Here we summarize the probe data by gene. For each gene, we use the mean M value of all probes annotated to that gene.
mvals.per.gene <- function(GRset, probe.gene.map, use.function=mean){
filtered.map <- lapply(probe.gene.map,
function(x){
intersect(x, rownames(GRset))
})
mvals <- getM(GRset)
out.list <- lapply(names(filtered.map),
function(ensg){
probe.data <- mvals[filtered.map[[ensg]], ]
if(length(nrow(probe.data)) ==0){
return(NULL)
}
return( apply(probe.data, 2, use.function, na.rm=TRUE) )
})
names(out.list) <- names(filtered.map)
out.tab <- do.call(rbind, out.list)
return(out.tab)
}
# This takes around 1 hour for each tissue, so it's commented out at the moment.
# mvals.per.gene.adipose <- mvals.per.gene(minfi.norm.a, probe.gene.map, use.function=mean)
# mvals.per.gene.muscle <- mvals.per.gene(minfi.norm.m, probe.gene.map, use.function=mean)
#
# write.csv(mvals.per.gene.adipose, file="gene_mean_mvals_adipose.csv", quote=FALSE)
# write.csv(mvals.per.gene.muscle, file="gene_mean_mvals_muscle.csv", quote=FALSE)
Finally, we export the M values as csv files. Since all probes overlppaoing SNPs have been filtered out, these data can be shared freely.
write.m.vals(minfi.norm.a, probes.near.genes, "./gene_mvals_adipose.csv.gz")
write.m.vals(minfi.norm.m, probes.near.genes, "./gene_mvals_muscle.csv.gz")
write.m.vals(minfi.norm.a, rownames(getM(minfi.norm.a)), "./non_snp_mvals_adipose.csv.gz")
write.m.vals(minfi.norm.m, rownames(getM(minfi.norm.m)), "./non_snp_mvals_muscle.csv.gz")
R.Version()
## $platform
## [1] "x86_64-apple-darwin15.6.0"
##
## $arch
## [1] "x86_64"
##
## $os
## [1] "darwin15.6.0"
##
## $system
## [1] "x86_64, darwin15.6.0"
##
## $status
## [1] ""
##
## $major
## [1] "3"
##
## $minor
## [1] "4.3"
##
## $year
## [1] "2017"
##
## $month
## [1] "11"
##
## $day
## [1] "30"
##
## $`svn rev`
## [1] "73796"
##
## $language
## [1] "R"
##
## $version.string
## [1] "R version 3.4.3 (2017-11-30)"
##
## $nickname
## [1] "Kite-Eating Tree"
sessionInfo(package = NULL)
## R version 3.4.3 (2017-11-30)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: OS X El Capitan 10.11.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] grid splines stats4 parallel stats graphics grDevices
## [8] utils datasets methods base
##
## other attached packages:
## [1] bindrcpp_0.2.2
## [2] gridExtra_2.3
## [3] pheatmap_1.0.8
## [4] biomaRt_2.34.2
## [5] EnsDb.Hsapiens.v75_2.99.0
## [6] ensembldb_2.2.2
## [7] AnnotationFilter_1.2.0
## [8] GenomicFeatures_1.30.3
## [9] AnnotationDbi_1.40.0
## [10] ChIPpeakAnno_3.12.7
## [11] VennDiagram_1.6.20
## [12] futile.logger_1.4.3
## [13] RColorBrewer_1.1-2
## [14] DMRcate_1.14.0
## [15] DMRcatedata_1.14.0
## [16] DSS_2.26.0
## [17] bsseq_1.14.0
## [18] limma_3.34.9
## [19] doParallel_1.0.11
## [20] sva_3.26.0
## [21] BiocParallel_1.12.0
## [22] genefilter_1.60.0
## [23] mgcv_1.8-23
## [24] nlme_3.1-137
## [25] minfiData_0.24.0
## [26] IlluminaHumanMethylation450kanno.ilmn12.hg19_0.6.0
## [27] IlluminaHumanMethylation450kmanifest_0.4.0
## [28] minfi_1.24.0
## [29] bumphunter_1.20.0
## [30] locfit_1.5-9.1
## [31] iterators_1.0.9
## [32] foreach_1.4.4
## [33] Biostrings_2.46.0
## [34] XVector_0.18.0
## [35] SummarizedExperiment_1.8.1
## [36] DelayedArray_0.4.1
## [37] matrixStats_0.53.1
## [38] Biobase_2.38.0
## [39] GenomicRanges_1.30.3
## [40] GenomeInfoDb_1.14.0
## [41] IRanges_2.12.0
## [42] S4Vectors_0.16.0
## [43] BiocGenerics_0.24.0
## [44] forcats_0.3.0
## [45] stringr_1.3.1
## [46] dplyr_0.7.5
## [47] purrr_0.2.4
## [48] readr_1.1.1
## [49] tidyr_0.8.1
## [50] tibble_1.4.2
## [51] ggplot2_2.2.1
## [52] tidyverse_1.2.1
## [53] data.table_1.10.4-3
## [54] knitr_1.20
##
## loaded via a namespace (and not attached):
## [1] R.utils_2.6.0
## [2] tidyselect_0.2.4
## [3] IlluminaHumanMethylationEPICanno.ilm10b2.hg19_0.6.0
## [4] RSQLite_2.1.1
## [5] htmlwidgets_1.2
## [6] munsell_0.4.3
## [7] codetools_0.2-15
## [8] preprocessCore_1.40.0
## [9] statmod_1.4.30
## [10] colorspace_1.3-2
## [11] BiocInstaller_1.28.0
## [12] rstudioapi_0.7
## [13] labeling_0.3
## [14] GenomeInfoDbData_1.0.0
## [15] mnormt_1.5-5
## [16] bit64_0.9-7
## [17] rprojroot_1.3-2
## [18] lambda.r_1.2.2
## [19] biovizBase_1.26.0
## [20] regioneR_1.10.0
## [21] R6_2.2.2
## [22] illuminaio_0.20.0
## [23] idr_1.2
## [24] bitops_1.0-6
## [25] reshape_0.8.7
## [26] assertthat_0.2.0
## [27] promises_1.0.1
## [28] scales_0.5.0
## [29] nnet_7.3-12
## [30] gtable_0.2.0
## [31] methylumi_2.24.1
## [32] rlang_0.2.0
## [33] rtracklayer_1.38.3
## [34] lazyeval_0.2.1
## [35] acepack_1.4.1
## [36] GEOquery_2.46.15
## [37] dichromat_2.0-0
## [38] broom_0.4.4
## [39] checkmate_1.8.5
## [40] yaml_2.1.19
## [41] reshape2_1.4.3
## [42] modelr_0.1.2
## [43] backports_1.1.2
## [44] httpuv_1.4.3
## [45] Hmisc_4.1-1
## [46] RBGL_1.54.0
## [47] RMySQL_0.10.15
## [48] tools_3.4.3
## [49] psych_1.8.4
## [50] nor1mix_1.2-3
## [51] siggenes_1.52.0
## [52] Rcpp_0.12.16
## [53] plyr_1.8.4
## [54] base64enc_0.1-3
## [55] progress_1.1.2
## [56] zlibbioc_1.24.0
## [57] BiasedUrn_1.07
## [58] RCurl_1.95-4.10
## [59] prettyunits_1.0.2
## [60] rpart_4.1-13
## [61] openssl_1.0.1
## [62] haven_1.1.1
## [63] cluster_2.0.7-1
## [64] magrittr_1.5
## [65] futile.options_1.0.1
## [66] ProtGenerics_1.10.0
## [67] missMethyl_1.12.0
## [68] hms_0.4.2
## [69] mime_0.5
## [70] evaluate_0.10.1
## [71] xtable_1.8-2
## [72] XML_3.98-1.11
## [73] mclust_5.4
## [74] readxl_1.1.0
## [75] compiler_3.4.3
## [76] crayon_1.3.4
## [77] R.oo_1.22.0
## [78] htmltools_0.3.6
## [79] later_0.7.2
## [80] Formula_1.2-3
## [81] lubridate_1.7.4
## [82] DBI_1.0.0
## [83] formatR_1.5
## [84] MASS_7.3-50
## [85] ade4_1.7-11
## [86] Matrix_1.2-14
## [87] permute_0.9-4
## [88] cli_1.0.0
## [89] quadprog_1.5-5
## [90] R.methodsS3_1.7.1
## [91] Gviz_1.22.3
## [92] bindr_0.1.1
## [93] pkgconfig_2.0.1
## [94] GenomicAlignments_1.14.2
## [95] registry_0.5
## [96] foreign_0.8-70
## [97] xml2_1.2.0
## [98] annotate_1.56.2
## [99] rngtools_1.3.1
## [100] ruv_0.9.7
## [101] pkgmaker_0.22
## [102] multtest_2.34.0
## [103] beanplot_1.2
## [104] rvest_0.3.2
## [105] doRNG_1.6.6
## [106] VariantAnnotation_1.24.5
## [107] digest_0.6.15
## [108] graph_1.56.0
## [109] rmarkdown_1.9
## [110] base64_2.0
## [111] cellranger_1.1.0
## [112] htmlTable_1.11.2
## [113] curl_3.2
## [114] shiny_1.0.5
## [115] Rsamtools_1.30.0
## [116] gtools_3.5.0
## [117] jsonlite_1.5
## [118] seqinr_3.4-5
## [119] BSgenome_1.46.0
## [120] pillar_1.2.2
## [121] lattice_0.20-35
## [122] GO.db_3.5.0
## [123] httr_1.3.1
## [124] survival_2.42-3
## [125] interactiveDisplayBase_1.16.0
## [126] glue_1.2.0
## [127] bit_1.1-13
## [128] stringi_1.2.2
## [129] blob_1.1.1
## [130] org.Hs.eg.db_3.5.0
## [131] AnnotationHub_2.10.1
## [132] latticeExtra_0.6-28
## [133] memoise_1.1.0